# Make rasters, TINs and Shapefiles from Lidar LAS files
#
# Author: Chris Harding, charding@iastate.edu
# Created: 2009 (cleaned up for the GIS workshop at MSU, June 2010)
# Last changed: Nov. 5, 2010 for Scott
# Copyright: use as you see fit, but some credit would be nice :)
#
#
# within a selected folder, converts all .las LIDAR data files into
# raster files and optional TIN and multipoint shapefiles.
# TINs can be decimated and their points stored in another shape file.
# The raster resolution and a elevation tolerance for the decimation are set
# by the user.
#
# This script can either be run within the folder that contains the las files to be
# converted or launched via a Tool from a ArcTools toolbox. For running it directly as
# a python script, you should copy the script into the folder (say c:\las_files\),
# open a shell (Windows PowerShell or cmd.exe), cd into the folder (cd c:\las_files)
# and, assuming python is installed to run from the command line, type:
#
# python Lidar_conversion_script.py <<< maybe some arguments >>>
#
# Arguments:
# below is a list of arguments for the script, note that they are actually args from
# ArcGIS, and should be strings even for numbers ("10 meters" not 10) and use lowercase
# "true" "false" strings (as opposed to Python's True/False booleans). The coordinate system is
# usually a .prj file but can also be a very longs string with the content of the file. The way the script
# works is to first init all parameters with a default (which you can manually change in this
# script) and to the overwrite them with any given argument values. So, if no args are given to
# the script it will use the default parameters given in this script in the Set default parameters
# section, which is handy if you want differnt parameters for different folders - just copy the
# script into the folder, edit it with the parameters you want for these las files and run with no args:
# python Lidar_conversion_script.py
# Of course, you can also give it some args but keep in mind that any arg not given will get
# the in script default value and you can't "skip" args. So if you only use the first two args (index 1 and 2):
# python Lidar_conversion_script.py "C:\\temp" "20 meters"
# the rest of the parameters (3 to 7) will default to whatever is hardcoded in the python script
# (Could be improved - I didn't see much point in doing a better job for arg parsing with getopts, etc.)
# 
# Input arguments:
# (argv[0] = name of python script, this is always given automatically already ...) 
# argv[1] = Folder - folder name, e.g. "C:\\temp" (or "." for "the folder the script is in")
#   of the folder containing the las. files you want to convert (e.g. l1.las, l2.las, etc.). All output
#   files/folders will also appear in this folder using the pre . part as "basename", e.g.
#   l1.shp for the las file's points as a Shapefile. Names of files/folders will use _ not spaces!
# argv[2] = Raster resolution (e.g. "10 meters"), in ERDAS imagine (float) format,
#   file name is las file name + DEM + resolution + .img, e.g. l1_DEM_10_Meters.img
# argv[3] = Decimation Z tolerance (e.g. "0.2 meters")  implies that the decimated TIN is saved as
#   las name + pts + decimation tolerance e.g. l1_tin_0.2_Meters . "0" is taken to mean no decimation  
# argv[4] = Save Decimated TIN's points in Shapefile? "true"/"false" (ignored if no decimation is set with 0!)
#   Shapefile name is las name + pts + decimation tolerance + .shp (etc.) e.g. l1_pts_0.2_Meters.shp
# argv[5] = Save Full TIN? "true"/"false"
#   TIN folder name is las name + tin e.g. l1_tin 
# argv[6] = Save Full TIN's points as Shapefile? "true"/"false"
#   Shapefile name is las name + pts + .shp e.g. l1_pts.shp
# argv[7] = Output Coordinate System
#   (e.g. "Coordinate Systems\Projected Coordinate Systems\Utm\Nad 1983\NAD 1983 UTM Zone 15N.prj" or
#   its long version:  "PROJCS['NAD_1983_UTM_Zone_15N',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-93.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]"

# function to choose where print outs go: either to AddMessage() or to stdout/stderr (or both)
# beware that it takes a sting, so you need to format your message before (I use %)
# Even if you're running this as a straight python script from a terminal, AddMessage should 
# print to the terminal, but Idle will need the python print instead.
def prnt(str):
    gp.AddMessage(str) # if you run this script from a shell you may want to comment out AddMessages and use print instead
    #print str
    #print >> sys.stderr, str # print to stderr instead of stdout
    
# ---------------------------------------------------------------------------
# Import system modules
import os, glob, sys
import platform # for checking which python version we're running

python_version_string = platform.python_version()
print "Running Python version", python_version_string

# Python 2.6 uses a different way to get the geoprocessing object
# than 2.5 (default for ArcGIS 9.3). Bail out for any other version of Python
if python_version_string[0:3] == "2.5":
    import arcgisscripting
    gp = arcgisscripting.create(9.3) # ArcGIS 9.3 geoprocessing 
    #gp = arcgisscripting.create() # pre 9.3 version
elif python_version_string[0:3] == "2.6":
    try:
        import arcgisscripting # if you use python 2.6, are you on ArcGIS 9.3 or 10?
        gp = arcgisscripting.create(9.3) #use the 9.3 arcscripting (still works for 10)

    except:
        # if not, maybe your're on ArcGIS 9.3 but instead of python 2.5 you use python 2.6 (naughty, naughty ...)?
        # install pywin32, see http://sourceforge.net/projects/pywin32/files/ and use this:
        try:
            import win32com.client
            gp = win32com.client.Dispatch("esriGeoprocessing.GpDispatch.1")
        except:
            prnt ("Can't create geoprocessing object via python 2.6 in ArcGIS 9.3")
            sys.exit()
else:
    prnt ("Can't create geoprocessing object - bailing out!")
    sys.exit()

# XML log file in C:\Users\<user name>\AppData\Roaming\ESRI\ArcToolbox\History (Vista)
# or C:\Documents and Settings\<user name>\Application Data\ESRI\ArcToolbox\History (XP)
gp.LogHistory = True  

# Check out any necessary licenses - you cannot run this script without a 3D Analyst license!
gp.CheckOutExtension("3D")

# Load required toolbox (different path on Vista? 64 bit?) - you may need to change that ...
# is this needed?
#gp.AddToolbox("C:/Program Files (x86)/ArcGIS/ArcToolbox/Toolboxes/3D Analyst Tools.tbx") # 64 bit?
#gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/3D Analyst Tools.tbx") # 32 bit

#
# Set default parameters 
#
folder = os.getcwd() # set to current location of the script, but could be a folder like "C:\\temp"
raster_resolution = "10 Meters" # cell size of raster
decimate_z_tolerance = "0.2 Meters"  #  or "0" for no decimation
save_decimated_points = "true" # "true" or "false"
save_full_TIN = "true"
save_full_points = "true"
coordinate_system = "Coordinate Systems\Projected Coordinate Systems\Utm\Nad 1983\NAD 1983 UTM Zone 15N.prj"

# if we have args, print them out and overwrite the defaults with them 
args = sys.argv
for i in range(1, len(args)):
    #prnt("arg %i is %s" % (i, args[i])) # debug, show args in ArcGIS tool "terminal"
   
    # not very elegant, should have used a dictionary for the parameters ...
    if i == 1: folder = args[1]
    if i == 2: raster_resolution = args[2]
    if i == 3: decimate_z_tolerance = args[3]
    if i == 4: save_decimated_points = args[4]
    if i == 5: save_full_TIN = args[5]
    if i == 6: save_full_points = args[6]
    if i == 7: coordinate_system = args[7]
 
gp.workspace = folder  # makes the folder with las files the official ArcGIS workspace
prnt("working folder is %s" % (folder)) 

# set coordinate system - I don't what the difference between 
# OutputCoordinateSystem and cartographicCoordinateSystem is, so I set both to be sure.
# However, I still had some problems with setting the CS, so I later set the CS again 
# directly on the Shapefile (not sure why).
gp.OutputCoordinateSystem = coordinate_system
gp.cartographicCoordinateSystem = coordinate_system

# allow output overwrite, unless this is done, it won't overwrite existing files!
gp.overwriteoutput = 1 # (Even when set it may not delete shapefiles)

# for all .las files in current folder do:
for LAS_file in glob.glob(folder + "\\*.las"):
    fname_prefix =  LAS_file[:-4]  #chop off the .las part of current file

    # make filenames for the various data files to eb created
    TIN_full_res = fname_prefix +  "_TIN"   # name of full resolution TIN (really a folder?)
    LIDAR_points_full_res = fname_prefix + "_pts.shp"  # full resolution points as multipoint shapefile 
    raster_resolution_unit = raster_resolution.replace(" ", "_") # replaces spaces with underscores
    DEM =  fname_prefix +  "_DEM_" + raster_resolution_unit + ".img"  # ERDAS imagine raster
    decimate_z_tolerance_unit = decimate_z_tolerance.replace(" ", "_") # replaces spaces with underscores
    TIN_dec = fname_prefix +  "_TIN_dec_" + decimate_z_tolerance_unit
    points_dec_shp = fname_prefix + "_pts_dec_" + decimate_z_tolerance_unit + ".shp"

    #
    # Process: LAS to Multipoint using LASToMultipoint_3d()
    #
    prnt("\nConverting %s to %s" % (LAS_file, LIDAR_points_full_res))
    
    # args to LASToMultipoint_3d are dependent on who wrote the LAS file, in my case:
    # - class code needs to be 2 (according to Chris Kahle from the Iowa DNR) for bare earth, which is different from last returns (?)
    # - return values should be ANY_RETURNS
    # - input_coordinate_system can be empty
    # - file_suffix is typically las
    # - I don't use attribute names
    # for more see ArcGIS Desktop Help for LAS to Multipoint (3D Analyst)
    average_point_spacing = "0.0000001" # to get all data points
    class_code = "2"  # or: "2;3", separate lists by ;
    return_values = "ANY_RETURNS" # or: "1;4;5"
    attribute_names = ""  # or: "INTENSITY Intensity;CLASSIFICATION Class;GPS_TIME GpsTime"
    input_coordinate_system = "" # means: take whatever is defined inside the las file
    file_suffix = "las"

    try:   
        res = gp.LASToMultipoint_3d(LAS_file, LIDAR_points_full_res, 
                        average_point_spacing, class_code, return_values, 
                        attribute_names, input_coordinate_system, file_suffix)
        gp.AddMessage(" result: " + res) # adds info about the processing to the file meta data (spatial - lineage)  
    except:
        prnt(gp.GetMessages())  # print error message
    
    # For me, the created shapefile sometimes still has no projection, 
    # so I explicitly set it here (again) with Define Projection:
    gp.toolbox = "management"
    gp.defineprojection(LIDAR_points_full_res, coordinate_system)
    
    #
    # Process: Create empty TIN and fill with points from shapefile
    #
    prnt("\nConverting shapefile %s to TIN %s" % (LIDAR_points_full_res, TIN_full_res))
    try:
        res = gp.CreateTin_3d(TIN_full_res, "")
        gp.AddMessage("CreateTin result: " + res)
    except:
        prnt(gp.GetMessages())
    
    try:
        res = gp.EditTin_3d(TIN_full_res, LIDAR_points_full_res + " Shape <None> masspoints true")
        # EditTin_3d takes only 2 args: the input TIN and another string with args e.g. "fullres.shp Shape <None> masspoints true"
        gp.AddMessage(" result: " + res)
    except:
        prnt(gp.GetMessages())
    
    # Decimate TIN ?
    tol = decimate_z_tolerance.split()[0] # get first part of string "0.1 Meters"
    if float(tol) > 0.0:  # a decimation tolerance of 0 means: don't decimate
        #
        # Process: Decimate TIN Nodes with decimate_z_tolerance
        #
        prnt("\nDecimating %s to %s" % (TIN_full_res, TIN_dec))
        try:
            res = gp.DecimateTinNodes_3d(TIN_full_res, TIN_dec, "ZTOLERANCE " + decimate_z_tolerance + " 0, \"NO_BREAKLINES\" ")
            # Again, 1 string for all args, note the \" in "NO_BREAKLINES"
            gp.AddMessage(" result: " + res) # doesn't create meta data?
        except:
            prnt(gp.GetMessages())
            
        #
        # Save decimated TIN nodes as point shapefile?
        #
        if save_decimated_points == "true":
            prnt("\nSaving decimated TIN %s as shapefile %s" % (TIN_dec, points_dec_shp))
            try:
                res = gp.TinNode_3d(TIN_dec, points_dec_shp, "", "Tag_Value")
                gp.AddMessage(" result: " + res)
            except:
                prnt(gp.GetMessages())    
                      
    #
    # Process: convert (full res) TIN to float ERDAS .img raster using 
    # natural neightbors interpolation.
    # z-factor of 1 assumes vertical unit is the same as horizontal
    #
    prnt("\nConverting full res TIN %s to raster %s" % (TIN_full_res, DEM))
    try:
        cellsize = "CELLSIZE " + raster_resolution
        z_factor = "1"
        res = gp.TinRaster_3d(TIN_full_res, DEM, "FLOAT", "NATURAL_NEIGHBORS", cellsize, z_factor)
        gp.AddMessage(" result: " + res)
    except:
        prnt(gp.GetMessages())

    # keep full res TIN around?
    if save_full_TIN == "false": # no, delete folder  
        gp.delete_management(TIN_full_res)

    # keep full resolution shapefile around? 
    if save_full_points == "false": # no, delete shapfile  
        gp.delete_management(LIDAR_points_full_res)
   
prnt("\nLidar conversion done")
#this prevents the cmd terminal to diappear after double clicking on the .py file but Arc will complain
#raw_input("done")





